First order phase transition in the Quantum Adiabatic Algorithm 
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We simulate the quantum adiabatic algorithm (QAA) for the exact cover problem for sizes up 
to N = 256 using quantum Monte Carlo simulations incorporating parallel tempering. At large N 
we find that some instances have a discontinuous (first order) quantum phase transition during the 
evolution of the QAA. This fraction increases with increasing N and may tend to 1 for N — s> oo. 
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It is of great interest to know if an eventual quantum 
computer could solve a broad range of hard "optimiza- 
tion" problems more efficiently than a classical computer. 
An important class is the NP-hard category [1] , for which 
it is believed that all classical algorithms take a time 
which grows exponentially with the problem size N. 

The most promising approach to solving optimization 
problems on a quantum computer seems to be the quan- 
tum adiabatic algorithm (QAA), proposed by Farhi et 
al. [2L The idea, which is related to quantum anneal- 
ing [3[, is that one adds to a "problem" Hamiltonian, 
Hp, whose ground state represents a solution of a clas- 
sical optimization problem, a non-commuting "driver" 
Hamiltonian, Hp>, so the total Hamiltonian is 



H{s) = (1 - s)H D + sHi 



(1) 



where s(t) is a time dependent control parameter, rip is 
expressed in terms of classical Ising spins taking values 
±1, or equivalently in terms of the z-components of the 
Pauli matrices for each spin, af. The simplest driver 
Hamiltonian is then then Hp, = — Yli=i <*t where of is 
the x-component Pauli matrix. 

The control parameter s(t) is at t = 0, so H—Hp,, 
which has a trivial ground state in which all 2 N basis 
states (in the a z basis) have equal amplitude. It then 
increases with t, reaching 1 at t = T, where T is the 
runtime of the algorithm, at which point H=Hp. If the 
time evolution of s(t) is sufficiently slow, the process will 
be adiabatic. Hence, starting the system in the ground 
state of Hp,, the system will end up in the ground state 
of Hp and the problem is solved. The time T required to 
find the ground state with significant probability is called 
the complexity. The bottleneck of the QAA is likely to 
be at one or more points where the energy gap from the 
ground state to the first excited state becomes very small, 
possibly due to a quantum phase transition. 

Early numerical work 0, 0] on very small systems, 
N < 24 (for a particular constraint satisfaction prob- 
lem known as "exact cover 3", also called l-in-3 sat), 
found that the complexity scaled polynomially with size, 
roughly as N 2 , which caused a good deal of excitement. 



However, this power law complexity may be an artifact 
of the very small sizes studied, so it is of great interest to 
determine whether the complexity continues to be poly- 
nomial for much larger sizes or whether a "crossover" to 
exponential complexity is seen. 

In previous work [5| we have used quantum Monte 
Carlo (QMC) simulations to investigate much larger sizes 
of the exact cover problem, up to N = 128. We found 
evidence that, while the median complexity is still poly- 
nomial, an increasing fraction of instances became very 
hard to equilibrate for the larger sizes. We have now con- 
siderably improved the algorithm, borrowing techniques 
from the the spin glass field to speed up equilibration. We 
have therefore been able to understand much better these 
"troublesome" instances, and find that they have a first 
order quantum phase transition. Furthermore, we have 
increased the range of sizes still further, up to N = 256, 
finding that the fraction instances with a first order tran- 
sition continues to increase with N, plausibly tending 1 
for N — > oo. The gap at a first order phase transition is 
likely to be exponentially small 0, @] , and hence lead to 
exponential complexity for the QAA. 

We now describe the model and our results in more 
detail. To make a comparison with the earlier work we 
study (essentially) the same model for Hp used by Farhi 
et al. [2]. This problem, known as exact cover, is a ran- 
dom satisfiability problem, a class which is known to be 
NP-hard. In exact cover there are N Ising spins and M 
"clauses" each of which involves three spins (chosen at 
random). The energy of a clause is zero if one spin is 
— 1 and the other two are 1, and otherwise the energy is 
a positive integer. The simplest Hamiltonian with this 
property is [8j 
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where ot\, cti and are the three spins in clause a, and 
the {of}(i = 1, ••• ,N) are Pauli matrices. In the ab- 
sence of the driver Hamiltonian, the Pauli matrices can 
be replaced by classical Ising spins Si taking values ±1. 
An instance has a "satisfying assignment" if there is at 
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N 


16 


32 


64 


128 


192 


256 


M 


12 


23 


44 


86 


126 


166 


Q 


0.7500 


0.7188 


0.6875 


0.6719 


0.6563 


0.6484 



We have considerably improved the algorithm relative 
to that in Ref. Q by incorporating "parallel temper- 



TABLE I: The sizes studied in the simulation. 

least one choice for the spins where the total energy is 
zero. As the ratio a = M/N is increased, there is a 
phase transition at a s where the number of satisfying 
assignments goes to zero. The version used by Farhi et 
al. considers only instances with a unique satisfying as- 
signment (USA), i.e. there is only one state with energy 
0. This has the advantage that the gap AE(s) between 
the ground state and first excited state is greater than 
zero in both limiting cases, H. = ~Hb and % = Hp, but 
will have a minimum at an intermediate value s = s* . In 
addition, it ensures that we work close to the satisfiabil- 
ity transition where the problem is particularly hard Q . 
Hence here, and in the earlier work Q, we consider in- 
stances with a USA. 

The method of generating instances with a USA is de- 
scribed in Ref. 5( . For each size N we choose the number 
of clauses M which maximizes the probability of finding a 
USA, see Table HI The actual number of spins simulated 
N', is somewhat less than N due to isolated sites being 
omitted, and others that do not affect the complexity are 
also "pruned off" The value of a = M/N seems to 
be close to the critical value a s ~ 0.626 [l(| for N — > oo. 

In QMC we simulate an effective classical model with 
Ising spins Si(r) = ±1 in which r (0 < r < j3 = T _1 ) 
is imaginary time. Following common practice, we dis- 
cretize imaginary time into L T "time slices" each repre- 
senting At = f3/L T of imaginary time. We take At = 1. 

As discussed previously @, for /3AE ^> 1 (where 
AE = Ei — E is the energy gap), and r <C (3, the 
time-dependent correlation function 

-, N' L T 

( 3 ) 

1=1 T = l 



C(r) 



N'L T 



is a sum of exponentials, i.e. 

C(T) = q + J2 A neM-(K-E )T}, (4) 



n>l 



where the A n are constants and q, the spin glass order 
parameter, is given by 



N< 



(5) 



At large t, the sum in Eq. (|4| is dominated by the term 
corresponding to the first excited state, {n = 1), and so 
AE = E\ — E can be obtained by fitting log[C(T) — q] 
against r for large t. 
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111 1 1 21 ] , which has been very successful in speeding 
up simulations of spin glass systems. Whereas in spin 
glasses, one simulates copies of the system at different, 
close-by temperatures, in the quantum case, the copies 
are at different values of the control parameter s. 

As already mentioned, the focus of the present study 
is to determine which instances have a first order tran- 
sition. Parallel tempering is very good at equilibrating 
the system on either side of the transition. However, it 
is still difficult (i) to determine exactly where the tran- 
sition occurs, because both phases are metastable in the 
region where they are not the equilibrium state, and (ii) 
to accurately determine the minimum gap for first-order 
instances, because it is so small. We have performed runs 
starting the spins both from a random initial configura- 
tion, and from the solution of the problem Hamiltonian. 
If we start by "seeding" the spins with the exact solu- 
tion, we are confident that the Monte Carlo is in the 
correct phase for s close to 1. It is also in the correct 
phase for small s because equilibration is easy in this re- 
gion. Hence, if a long simulation starting the spins from 
the exact solution produces a sharp discontinuity, we feel 
that this is almost certainly the correct behavior. 

In order to investigate whether or not a first order tran- 
sition occurs we focus on the spin glass order parameter 
q defined in Eq. ([5]). The expectation value of q is al- 
ways non-zero because of terms linear in the a z (magnetic 
field terms) in the Hamiltonian, Eq. ([2]). To determine 
the square of the average without bias we simulated two 
copies of the spins at each value of s and evaluate (Si) 2 
as (Si)^ {Si)( 2 K A representative result for an instance 
with a first order transition is shown in Fig. [TJ Note 
the very rapid increase in q over a very small range of 
s, and that the two curves on each side of the jump are 
obviously displaced vertically with respect to each other. 

The dip before the jump clearly seen in Fig. Q] provides 
clear evidence for a two-phase coexistence, and hence a 
first order transition, for the following reason. If both 
copies are in the same phase, then the mean value of 
(Si) is the same in both copies. However, right at the 
first order transition, one copy can be in one phase (the 
low-<7 phase, say) and the other copy in the other (high- 
q) phase. The average value of (Si) can have different 
signs in the two phases for some sites i. Hence, the typ- 
ical Hamming distance between the spin configurations 
in the two copies can be even greater (and so q even 
smaller) than when both copies are in the low-g phase. 
In every instance where we observed a sharp jump, this 
was proceeded by a dip. Hence we use the dip as a precise 
criterion for a first order transition. 

Of course, even a first order transition is rounded out 
for a finite-size system. To estimate the size of the round- 
ing we need to consider the two cases AE m i n ^> T and 
AE min <C T separately, where AE mnl is the minimum 
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FIG. 1: (Color online) An instance with a first order transition 
with N = 128. Note the expanded horizontal scale. 

value of the gap at the transition. If AE min ^> T, 5s 
is the range of s over which AE changes by an amount 
AE mm , whereas if AE mul <C T, Ss is the range of s over 
which AE changes by an amount equal to T. Hence 

f AE min (^fy 1 , (AE min » T), 
Ss = < (6) 

l T (^)~' (A£ mi „«T). 

Figure [2] shows the finite-size rounding for an instance 
with N = 64, small enough that we can equilibrate 
through the (first order) transition. For j3 < 1024 the 
width of the transition region increases as /3 = 1/T de- 
creases, but for (3 > 1024 the width is independent of j3. 
For this instance we find AE nl [ n — 0.0021 as shown in the 
inset, so the width of the rounding becomes independent 
of T when T <C AE as expected. 

In Fig. |3] we plot the fraction of instances with a first 
order transition. For each size we have studied N- lnst — 50 
instances. If we denote the first-order fraction by r then 
the error bar in r is y/r{l — r)/(A r ; nst — 1). The figure 
shows that r increases rapidly with N and, very plausi- 
bly, tends to 1 for N —> oo. We see that the first order 
fraction is slightly greater than a half for N = 128. In our 
earlier work [5j we found that the median complexity con- 
tinued to be polynomial up to N — 128 (the largest size 
studied). However, there is no contrast with the present 
work because, as already noted in Ref. Q, the models 
used are slightly different, and as a result the crossover 
to a first order transition occurs at a slightly lower value 
of N in the present model. The crossover to first order 
would have been seen in the earlier model if somewhat 



FIG. 2: (Color online) The main figure shows the spin glass 
order parameter q, defined in Eq. as a function of s for an 
instance with TV = 64 which has a first order transition. The 
different curves are for different values of /3. The inset shows 
the energy gap AE as a function of s for f3 — 2048, indicating 
that A75 m i n = 0.0021 (same value was found for /3 = 1024 
and 4096). From the main figure one sees that the width of 
the finite-size rounding increases with T = 1//9 for T 3> AE 
but is independent of T in the opposite limit T <C AE, as 
expected from Eq. ([6]). Note the expanded horizontal scale. 
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FIG. 3: (Color online) The fraction of instances with a first 
order transition (defined in the way discussed in the text) as 
a function of size. For each size, 50 instances were studied. 
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larger sizes had been studied. 

Exponentially small gaps have been discussed before 
in the context of the QAA. Some time ago, one of us [l3| 
pointed out for a different problem, number partitioning, 
that the minimum gap is exponentially small, because of 
a transition between the states that are "localized" and 
"extended" in the computational basis. 

Altshuler et al. Q predict an exponentially small gap 
at large N for exact cover. Performing perturbation the- 
ory away from s = 1 they argue that there will be a level 
crossing between two "localized" states for s close to 1 
at which point the ground state configuration changes 
abruptly. In our numerics, there is a big variation in the 
location of the first order transition for a given size, but 
we do not detect a systematic shift towards s = 1 as 
the size increases. However, Altshuler et al. predict that 
1 — s ~ N^ 1 / 8 , which is probably too slow to be visible 
in our data. It will be interesting to investigate in future 
work whether the first order transition found here is due 
to the mechanism they propose. 

Farhi et al. [3] used a continuous imaginary time QMC 
method to study a very similar problem to ours, ex- 
cept that two solutions far away in Hamming space are 
"planted" into the Hamiltonian. This ensures that there 
is a finite probability of a first order transition where the 
equilibrium state changes from one planted solution to 
another. By contrast, our work does nothing explicit to 
impose a first order transition. 

Jorg et al. [HI studied quantum annealing for the quan- 
tum random energy model (REM) , the classical version of 
which has a "1-step replica symmetry breaking" (also 
called a "random first order" ) transition. Following Gold- 
schmidt [T3|, they find a discontinuous quantum transi- 
tion and argue that this leads to an exponentially small 
gap. They also observed that an exponentially small gap 
is seen in quantum versions of several models with ran- 
dom first order transitions and suggested that this may 
be the general feature of all such models, including sat- 
isfiability III 18| . However, the classical REM has zero 
spin glass order parameter q in the disordered phase [l6| 
whereas classical random satisfiability models have lower 
symmetry because q is always non-zero due to the terms 
linear in a z in Eq. ([2} . Consequently, it is not obvious to 
us that the first order quantum transition observed here 
is due to the same mechanism as that found [H| for the 
quantum REM. Very recently, a first order transition has 
also been found in another model by Jorg et al. . 

To conclude, we have a found a crossover to a first 
order quantum phase transition during the evolution of 
the QAA for instances of exact cover with a unique sat- 
isfying assignment when the size becomes greater than 
about 100. It is possible that the complexity for ran- 
dom instances of exact cover could be different. We are 
therefore studying instances of exact cover with the USA 
constraint removed, and will also study other models in 
addition to exact cover. 
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